%% 机器学习与数据预测——信号处理
% Machine Learning and Data Prediction -- Signal Processing
%% 1.泰勒级数应用Taylor series application
%%  1.1.diff,int,cumsum函数用法
% 符号函数求导(偏导),求积分(偏积分)
syms x y z m n
f=x^3+y^2+z^2+m^2+n^2;
y1=diff(f,x,1);y2=diff(f,y,2);y3=diff(f,z,3);
y4=int(f,x,1);y5=int(f,y,2);y6=int(f,z,3);
% 匿名函数求导
syms x1 x2 x3
f=@(x1,x2,x3)x1*x2*x3+x2^2+x3^2;
y1=diff(f,x1,1);y2=int(f,x2,2);
% 数组求差分
x1=[1,2,3,4,5,6];y1=diff(x1,1);y2=cumsum(x1);%此时用累加和代替int
% 矩阵：对列向量求差分
A=[1,2,3;4,5,6;7,8,9];y1=diff(A,1);y2=cumsum(A);
%%  1.2.对带有边界的函数求其泰勒展开式
% 符号函数用法：
% taylor(f,x,a,'order',n); %求f关于自变量x在x=a处的泰勒展开式前n项，若a缺省，则默认a=0
syms x1 x2 x3
f=x1^2+x2*x3+x1*x3;
y1=taylor(f,x1,1,'order',4);
y2=taylor(f,[x1,x2,x3],[1,2,1],'order',3);
% 一元数值函数编写：[A,coefficient]=Taylor_series(f,x0,n,xmin,dx,xmax)
% 函数说明：
% 输入量：f为函数信号，为数组；x0为泰勒展开的自变量取值点；n为泰勒级数阶数；
% xmin为自变量取值最小值；xmax为自变量取值最大值，dx为取值间隔。
% 输出量：A为各阶泰勒级数，A第i行代表第i阶泰勒展开式；coefficient为各阶泰勒
% 级数系数，第i项代表第i阶泰勒展开式系数
% 函数注意：
% n值不能过大，5以上效果较优；
% 输入f量不能过大。
n=6;%泰勒级数阶数
x0=0;%求级数的点位置
xmin=-10;xmax=10;dx=1e-4;
x=xmin:dx:xmax;
f=exp(x);
[A,coefficient]=Taylor_series(f,x0,n,xmin,dx,xmax);
%% 2.傅里叶级数与傅里叶变换Fourier series and Fourier transform
%%  2.1.信号基础
%%   2.1.1.对信号实施的基础运算：
a=2;
xmin=-10;xmax=10;dx=1e-4;
x=xmin:dx:xmax;
f=sin(x);
y1=a*f;plot(x,y1);% 数乘
y2=cumsum(f*dx);plot(x,y2);% 积分(不带常数项)
y3=diff(f)/dx;plot(x,y3);% 微分(求导)
t0=1;
x1=x+t0;f=sin(x);plot(x1,f);% 信号右移t0
x0=pi/2;x2=-x+2*x0;f=sin(x);plot(x2,f);% 信号绕点x0翻转
a=3;x3=a*x;plot(x3,f);% 信号缩展，其中a>0,a>1时信号展开；a<1时信号收缩
% 对变换之后的信号进行整理，使自变量从小到大变化
% 对[信号绕点x0翻转]举例：
% 传统方法：
m0=[f;x2];%将函数和自变量放在同一矩阵中，第一行为函数值，第二行为自变量取值
x2=sort(x2);k=zeros(1,length(x2));
for i=1:length(x2)
    [~,k0]=find(m0==x2(i));
    k(i)=k0;
end
f=f(k);plot(x2,f);
% 利用sort函数改进方法：(在matlab中推荐使用，可节省时间)
[x2,k]=sort(x2);
f=f(k);plot(x2,f);
%%   2.1.2.信号关于x0点处奇偶分解：
xmin=-10;xmax=10;dx=1e-4;x=xmin:dx:xmax;f=exp(x);% 原始信号
x0=0;x2=-x+2*x0;[x2,k]=sort(x2);f0=f(k);% 结合上述方法获得翻转信号
xnew_min=max([x2(1),x(1)]);xnew_max=min([x2(end),x(end)]);% 获得自变量范围
xnew=xnew_min:dx:xnew_max;
% 自变量范围为两信号范围求交集
[~,k1]=find(x>=xnew_min&x<=xnew_max);f1=f(k1);% 获得新自变量范围下原始信号
[~,k2]=find(x2>=xnew_min&x2<=xnew_max);f2=f0(k2);% 获得新自变量范围下原始信号
c0=min([length(f1),length(f2),length(xnew)]);% 简单检验信号长度，防止信号长度
% 不同导致程序无法运行
xnew=xnew(1:c0);
Odd_signal=1/2*(f1(1:c0)-f2(1:c0));% 奇信号
Even_even_signal=1/2*(f1(1:c0)+f2(1:c0));% 偶信号
plot(xnew,Odd_signal);hold on;plot(xnew,Even_even_signal);hold off;
legend('Odd_signal','Even_even_signal');
%%  2.2.傅里叶级数
%%   2.2.1.基础信号特征：
xmin=-10;xmax=10;dx=1e-4;x=xmin:dx:xmax;f=exp(x);% 原始信号
E=sum(f.^2*dx);% 求信号能量
P=E/(xmax-xmin);% 求信号功率
% 连续情况下，能进行傅里叶级数，傅里叶变换的信号需满足：
% 信号有有限个不连续点；有有限个极值点，满足绝对可积
% 一般情况下，所有数值处理的有限信号必满足以上条件
% 采样定理：对连续信号进行采样的频率需大于连续信号最高频率的两倍
%%   2.2.2.三角形式傅里叶级数：
xmin=-10;xmax=10;dx=1e-4;x=xmin:dx:xmax;f=exp(x);% 原始信号
% 将信号分解为三角形式，即：f(x)=a(0)+sum(a(k)*cos(k*w0*t)+b(k)*sin(k*w0*t))
% 编写数值计算函数：function [ak,bk]=Tri_Fourier_series(f,w0,km,xmin,dx,xmax)
% 函数说明：
% 输入量f为信号，w0为f基频，km为求傅里叶级数最大基频倍数，xmin为自变量最小值，
% dx为自变量变化间隔，xmax为自变量最大值，输出[ak,bk]为傅里叶级数系数。
xmin=-10;xmax=10;dx=1e-4;w0=2*pi;x=xmin:dx:xmax;
f=sin(w0*x)+sin(3*w0*x)+5*sin(5*w0*x);% 原始信号
km=10;
[ak,bk]=Tri_Fourier_series(f,w0,km,xmin,dx,xmax);
E=sum(f.^2*dx);% 求信号能量
P=E/(xmax-xmin);% 求信号功率
% 各次谐波含量可表示为：
A0=ak(1)^2/P; % 常数偏置项
A=(bk.^2+ak(2:end).^2)/(2*P); % 各次波分量
A=[A0,A];
%%   2.2.3.傅里叶变换：
% 符号运算
w0=2*pi;syms x w;
f=1/(1/(2*pi/w0)*x+1);
F=fourier(f,x,w);% 对自变量为x的函数f进行傅里叶变换，变换后自变量为w
f1=ifourier(F,w,x);% 对自变量为w的函数F进行傅里叶逆变换，变换后自变量为x
A=abs(f1);fplot(A);
% 数值运算
% 编写函数方法：
% 信号傅里叶变换为：F=int(f*exp(-1i*w*x)*dx)
% 编写数值计算函数：[F,A]=Fourier_transform(f,xmin,dx,xmax,wmin,dw,wmax)
% 函数说明：
% 输入量f为信号，xmin为自变量最小值，dx为自变量变化间隔，xmax为自变量最大值，
% wmin为频域最小值，dw为频域间隔，wmax为频域最大值，输出[F,A]为傅里叶变换后
% 函数与其对应的函数增益
% 注意：利用该函数求解，计算量为：length(wmin:dw:wmax)*length(xmin:dx:xmax)
% 则需要注意计算量不能过大，由于函数为数值求解方法，对于有确定周期函数有可能
% 计算不到周期对应的频率，因此该函数可定性用于对某非信号的频域分析
xmin=-10;% 自变量最小值
xmax=10;% 自变量最大值
dx=1e-3;% 自变量变化间隔
x=xmin:dx:xmax;
w0=2*pi;
f=sin(w0*x)+sin(3*w0*x)+sin(5*w0*x);% 原始信号
wmin=0;% 频域最小值
wmax=10*w0;% 频域最大值
dw=0.1;% 频域变化范围
[F,A]=Fourier_transform(f,xmin,dx,xmax,wmin,dw,wmax);
plot(wmin:dw:wmax,A);
% matlab内置函数：fft(快速傅里叶变换)
xmin=-10;% 自变量最小值
xmax=10;% 自变量最大值
dx=1e-4;% 自变量变化间隔(可视为采样频率)
T=1/dx;% 自变量变换周期(可视为采样周期)
w0=2*pi;
x=xmin:dx:xmax;
L=length(x);% 自变量个数
f=sin(w0*x)+sin(3*w0*x)+5*sin(5*w0*x);% 原始信号
F=fft(f);% 计算信号傅里叶变换
A2=abs(F/L);% 双侧频谱
A1=A2(1:round(L/2)+1);
A1(2:end-1)=2*A1(2:end-1);% 得到单侧频谱
w=dx*(0:round(L/2))/L;
plot(w,A1);
%% 3.拉普拉斯变换与线性系统设计Laplace Transform and Linear system design
%%  3.1.拉普拉斯变换
% 符号运算
w0=2*pi;syms x s;
f=1/(1/(2*pi/w0)*x+1);
L=laplace(f,x,s);% 对自变量为x的函数f进行拉普拉斯变换，变换后自变量为s
f1=ilaplace(L,s,x);% 对自变量为s的函数F进行逆拉普拉斯变换，变换后自变量为x
% 数值运算
% 信号拉普拉斯变换可表示为：F=int(f*exp(-s*x)*dx)
% 确定待求解拉氏变换自变量范围s
pmin=0;dp=1e-2;pmax=5;%实部范围
qmin=0;dq=1e-2;qmax=5;%虚部范围
p=pmin:dp:pmax;q=qmin:dq:qmax;
s=zeros(length(q),length(p));
L=zeros(length(q),length(p));
for i=1:length(p)
    s(:,i)=1i*q'+p(i);
end
% 原始信号
xmin=-10;xmax=10;dx=1e-4;w0=2*pi;x=xmin:dx:xmax;
f=sin(w0*x)+sin(3*w0*x)+5*sin(5*w0*x);
% 求解拉氏变换函数
% Laplace_transform=@(f,s,x,dx)sum(f.*exp(-s.*x)*dx);
for i=1:length(p)
    tic
    for j=1:length(q)
        L(j,i)=sum(f.*exp(-s(j,i).*x)*dx);
    end
    toc
end
LA=abs(L);
%%  3.2.线性系统设计
%%   3.2.1.系统设计基础方法：
% 传递函数：某线性定常系统，在零初始条件下，输出拉氏变换与输入拉氏变换之比
% 设传递函数为：G(s)，有典型环节传递函数：
% 比例环节：G(s)=K；惯性环节：G(s)=K/(T*s1)；
% 振荡环节：G(s)=K/(T^2*s^2+2*zeta*T*s+1)；积分环节：G(s)=K/s；
% 理想微分环节：G(s)=K*s；延迟环节：G(s)=K*exp(-ta*s)
% 原始信号
xmin=0;xmax=10;dx=1e-4;w0=2*pi;x=xmin:dx:xmax;
f=sin(w0*x)+sin(3*w0*x)+5*sin(5*w0*x);
% 建立线性定常系统一般形式传递函数：
num=[4,5,6];% 分子多项式系数
den=[1,2,3,4,5];% 分母多项式系数
G1=tf(num,den);% 建立传递函数结构体
y=lsim(G1,f,x);% 绘制输出信号，对于实践问题，尽可能让开始时间=0
% 建立线性定常系统零极点形式传递函数：
kgain=10;% 增益
Z=[-1,-2,-3+2*1i];% 分子零点
P=[-3,-7,-4,-2+2*1i,-3-3*1i];% 分母极点
G2=zpk(Z,P,kgain);% 建立传递函数结构体
% 两种形式传递函数可互换
G3=zpk(G1);G4=tf(G2);
% 传递函数的组合：G1+G2，G1*G2
G5=tf(G1+G2);G6=tf(G1*G2);
% 反馈连接结构传递函数：设G1为前向通道传递函数，G2为反馈通道传递函数，则：
G7=feedback(G1,G2,-1);% -1表示负反馈，+1表示正反馈，返回形式默认为零极点形式
z=tzero(G1);p=poly(G1);% 求取传递函数G1零极点
% 利用matlab simulink可设计线性定常系统
% 例：Simple_linear_timeinvariant_system.slx
% 模型初始化：原始信号
xmin=0;xmax=10;dx=1e-4;w0=2*pi;x=xmin:dx:xmax;
f=sin(w0*x)+sin(3*w0*x)+5*sin(5*w0*x);
f1=[x',f'];% simin矩阵形式，其中第一列为时间戳，第二列为信号数据
sim('Simple_linear_timeinvariant_system.slx')% 执行仿真文件，输出结果将保存
% 至out.out中，为数组形式
%%   3.2.2.简单系统性能分析：
% 对于一个稳定的线性系统，其对于暂态响应或冲激响应，最终输出趋近于0.
% 对于稳定的线性系统，通常采用单位阶跃响应作为测试信号，测试以下动态性能指标：
% 设输入单位阶跃信号为:
xmin=0;xmax=20;dx=1e-4;x=xmin:dx:xmax;
f=[1,ones(1,length(x)-1)];
% 输出为y，传递函数：
num=1;den=[1,1,1];G1=tf(num,den);% 建立传递函数
y=lsim(G1,f,x);% 输出
plot(x,y);
% 利用matlab函数step,impulse可以方便地求出系统的单位阶跃响应和冲激响应
step1=step(num,den);impulse1=impulse(num,den);
% 超调量theta：theta=(max(y)-y(end))/y(end);反映系统相对稳定性，超调量越接近0，相对稳定性越好
theta=(max(y)-y(end))/y(end);
% 调节时间ts：ts满足abs((y(round((ts-xmin)/dx+1))-y(end))/y(end))<=delta，其中，
% delta=0.02或0.05，为规定的允许误差
delta=0.019;
for ts=sort(x,'descend')
    if(abs((y(round((ts-xmin)/dx+1))-y(end))/y(end))>delta)
        break
    end
end
% 延迟时间td：系统阶跃响应达到稳态值50%时间，即满足：y(round((ts-xmin)/dx+1))=0.5*y(end)
for td=x
    if(abs(y(round((td-xmin)/dx)+1)-0.5*y(end))<1e-4)
        break
    end
end
% 上升时间tr：系统阶跃响从稳态值10%到第一次达到稳态值90%所需时间，
% 即满足：y(round((tr1-xmin)/dx)+1)=0.1*y(end);y(round((tr2-xmin)/dx)+1)=0.9*y(end);
% tr=tr2-tr1;
for tr=x
    if(abs(y(round((tr-xmin)/dx)+1)-0.1*y(end))<1e-4)
        tr1=tr;
    end
    if(abs(y(round((tr-xmin)/dx)+1)-0.9*y(end))<1e-4)
        tr2=tr;
        break
    end
end
tr=tr2-tr1;
%%   3.2.3.简单系统根轨迹法和频率法分析：
% 对于具有负反馈的线性系统，其传递函数可表示为：W(s)=G(s)/(1+G(s)*H(s));
% 其中：G(s)为前向通道传递函数，H(s)为反馈通道传递函数，则有根轨迹基本方程：
% G(s)*H(s)=-1，称G(s)*H(s)为系统开环传递函数
% 已知连续系统开环传递函数为：G(s)*H(s)=K*(2*s^2+5*s+1)/(s^2+2*s+3);
num=[2,5,1];den=[1,2,3];
rl=rlocus(num,den);rlocus(num,den)% 绘制根轨迹图
sgrid;rlocus(num,den);
[K,P]=rlocfind(num,den);% 该命令可以在根轨迹图出现一光标，点击根轨迹图一点可以
% 返回该点极点坐标值P和对应增益K
% 对于某系统，当其输入为具有一定频率正弦信号时，其输出一般也为具有一定频率的
% 正弦信号，此时传递函数幅值相位会随频率发生变化，可以分别称之为幅频特性
% 和相频特性，其反映了系统对于不同频率信号的传递性能。
num=[2,5,1];den=[1,2,3];
wmin=0;wmax=1e3;dw=1;w=wmin:dw:wmax;% 频率范围
[mag,phase,w]=bode(num,den,w);% 返回传递函数在对应频率下幅值mag，相角phase
bode(num,den) % 绘制bode图
nyquist(num,den)% 绘制nyquist图
%% 4.Z变换与线性离散系统分析Z transform and linear discrete system analysis
%%  4.1.z变换
% 在以上的数值matlab分析中，从本质讲都是对连续系统离散化进行近似数值计算。
% 在本节讨论因果离散时间信号(自变量起始点为0)。
% 符号运算
syms n z;T=1e-3;% 采样周期
f=sin(T*n);
Z=ztrans(f,n,z);% 对自变量为n的序列f进行z变换，变换后自变量为z
f1=iztrans(Z,z,n);% 对自变量为z的函数Z进行逆z变换，变换后自变量为n
% 数值运算
% 信号z变换可表示为：F=sum(f.*z.^(-n))
% 确定待求解z变化自变量范围z
pmin=0;dp=1e-2;pmax=5;%实部范围
qmin=0;dq=1e-2;qmax=5;%虚部范围
p=pmin:dp:pmax;q=qmin:dq:qmax;
z=zeros(length(q),length(p));
Z=zeros(length(q),length(p));
for i=1:length(p)
    z(:,i)=1i*q'+p(i);
end
% 原始信号
T=1e-3;% 采样时间
tm=10;% 最大时间
t=T:T:tm;% 时间自变量范围
n=cumsum(ones(1,length(t)));
f=sin(t)+sin(3*t)+5*sin(5*t);
% 求解z变换
% z_transform=@(f,z,n)sum(f.*z.^(-n));
for i=1:length(p)
    tic
    for j=1:length(q)
        Z(j,i)=sum(f.*z(j,i).^(-n));
    end
    toc
end
ZA=abs(Z);
%%  4.2.线性离散系统设计
%%   4.2.1.离散系统设计基础方法：
% 原始信号
T=1e-3;% 采样时间
tm=10;% 最大时间
t=0:T:tm;% 时间自变量范围
n=cumsum(ones(1,length(t)));
f=sin(t)+sin(3*t)+5*sin(5*t);
% 建立线性离散系统一般形式传递函数：
num=[4,5,6];% 分子多项式系数
den=[1,2,3,4,5];% 分母多项式系数
G1=tf(num,den,T);% 建立离散系统传递函数结构体
y=lsim(G1,f,t);% 绘制输出信号，对于实践问题，尽可能让开始时间=0
% 建立线性离散系统零极点形式传递函数：
kgain=10;% 增益
Z=[-1,-2,-3+2*1i];% 分子零点
P=[-3,-7,-4,-2+2*1i,-3-3*1i];% 分母极点
G2=zpk(Z,P,kgain,T);% 建立离散传递函数结构体
% 反馈连接结构离散传递函数：设G1为前向通道离散传递函数，G2为反馈通道离散传递函数，则：
G=feedback(G1,G2,-1);% -1表示负反馈，+1表示正反馈，返回形式默认为零极点形式
% 连续系统传递函数变为离散系统传递函数
% 建立线性连续系统一般形式传递函数：
num=[4,5,6];% 分子多项式系数
den=[1,2,3,4,5];% 分母多项式系数
G1=tf(num,den);% 建立连续系统传递函数结构体
T=1e-3;% 采样时间
G3=c2d(G1,T);% 转化
% 离散系统传递函数变为连续系统传递函数
G4=d2c(G3);
% 利用matlab simulink可设计离散定常系统，其设计方法与连续系统相似，此处不再赘述
% 示例可见：Simple_discrete_linear_timeinvariant_system.slx
% 原始信号
T=1e-3;% 采样时间
tm=10;% 最大时间
t=0:T:tm;% 时间自变量范围
n=cumsum(ones(1,length(t)));
f=sin(t)+sin(3*t)+5*sin(5*t);
f1=[t',f'];
sim("Simple_discrete_linear_timeinvariant_system.slx")
% 对于一个稳定的线性离散系统，其对于暂态响应或冲激响应，最终输出趋近于0.
% 线性定常离散系统稳定的充要条件是系统闭环传递函数全部极点均分布在平面上以
% 原点为圆心的单位圆内，或者系统所有特征根的模均小于1。
%% 4.小波变换Wavelet transform
% 本节及之后经验模态分解将不再单独编写数值计算函数，针对不同问题，各部分内容
% 都有大量内容值的细究，限于篇幅，本文重点介绍其matlab基础相关函数使用方法。
% 本文意在给读者提供基础的matlab分析思路。
% 小波变换可以看作一种特殊的加窗傅里叶变换，传统傅里叶分析重点针对具有周期性
% 的信号，对于具有有限物理量的逻辑分析的数据具有广泛的应用价值，当某一系统影
% 响因素过多，影响大小及作用程度不确定时，系统输出信号往往表现出明显的非周期
% 性，趋势性，波动性和随机性等，对于该类信号的处理，小波变换提供了不错的处理
% 思路。
% 小波变换基本公式为：WT(a,tao)=1/sqrt(a)*int(f(t)*phi((t-tao)/a)*dt)
% 式中，函数phi为小波基函数，是具有一定时域和频率的函数，tao为位移因子，a为
% 尺度因子，对phi进行伸缩和位移变换，即可捕捉到信号不同时间处信号的频率变化。
%%  4.1.matlab简单小波变换基础
% 原始信号
T=1e-3;% 采样时间
tm=10;% 最大时间
t=0:T:tm;% 时间自变量范围
f=exp(-t/tm).*(sin(t)+sin(3*t)+sin(20*t));
wname1='morse';% 指定小波变换使用函数，可以为morse、Morlet (Gabor)和bump小波
wt=cwt(f,wname1,1/T);% 进行小波变换
cwt(f,wname1,1/T) % 绘制图像
f1=icwt(wt,wname1);% 进行小波逆变换
% 离散小波变换
wname2='db4';% 指定离散小波变换使用函数，可以为：
% morl(Morlet小波),mexh(墨西哥草帽小波),meyr(Meyer小波),haar(Haar小波),
% dbN(紧支集正交小波，其中N为正整数，之后相同),symN(近似对称紧支集正交小波),
% coifN(Coiflet小波),biorNr.Nd(双正交样条小波)
[cA,cD]=dwt(f,wname2);% 进行离散小波变换，返回近似信号cA(可对应低频)和细节信
% 号cD(可对应高频)
f2=idwt(cA,cD,wname2);% 进行逆离散小波变换
% 多层小波分解
N=3;% 3层小波分解
[c,l]=wavedec(f,N,wname2);% 进行N层小波分解
for i=1:N
    ap=appcoef(c,l,wname2,i);
    eval(['ap',num2str(i),'=ap;']);% 提取尺度为i(i=1,2,...,N)的低频系数
end
for i=1:N
    de=detcoef(c,l,wname2,i);
    eval(['de',num2str(i),'=de;']);% 提取尺度为i(i=1,2,...,N)的高频系数
end
f1=waverec(c,l,wname2);% 多层小波重构
%%  4.2.matlab小波降噪
% 初始信号，这里采用matlab内置信号noisdopp
load noisdopp;f=noisdopp;
wname3='db4';% 指定小波降噪使用函数，matlab常用内置有：
% bl(Daubechies),beyl(Beylkin),Coif(Coiflets),db(Daubechies),Haar(Haar),
% sym(Symlets),vaid(Vaidyanathan)等
method='Bayes';% 指定降噪方法，matlab常用内置有：
% Bayes,BlockJS,FDR,Minimax,SURE,UniversalThreshold
f1=wdenoise(f,Wavelet=wname3,DenoisingMethod=method);
%% 5.经验模态分解法及其拓展Empirical mode decomposition method and its extension
%%  5.1.经验模态分解
% 经验模态分解常应用在非平稳信号的局部频谱分析上，相比较傅里叶变换，小波变换
% 经验模态分解利用信号中极值点的信息对信号进行分解，其对于频率的捕捉通过信号
% 极值点信息进行反映，通过经验模态分解算法，可以将信号分解为若干本征模态函数
% （IMF）和具有单调性的残差（residual），此时，信号可表示为：
% f(t)=IMF1f(t)+IMF2f(t)+...+IMFnf(t)+residual(t);
% matlab函数如下：
% 初始信号，这里采用matlab内置信号
load('sinusoidalSignalExampleData.mat','X','fs')
t = (0:length(X)-1)/fs;% fs为采样频率，X为初始信号
[imf,residual]=emd(X);% 进行EMD分解
hht(imf,fs) % 绘制imf分量频谱图
%% 拓展：变分模态分解，集合经验模态分解
% 信号进行EMD分解后，同一频率的量可能出现在不同imf分量中，即出现模态混叠现象
% 同时EMD对信号边界点处处理具有局限性，变分模态分解，集合经验模态分解可以较好
% 的解决以上问题。
% matlab函数如下：
[Bimf,residualv]=vmd(X);% 进行VMD分解
% matlab中没有内置eemd分解函数，可以下载台湾中央大学emd工具箱进行编写，下载
% 地址为：http://khsci.com/docs/index.php/2020/04/09/1/
% 下边介绍该工具箱中eemd分解部分
Nstd=0.2;% Nstd为附加噪声标准差与Y标准差之比
NE=100;% NE为对信号的平均次数
imf1=eemd(X,Nstd,NE);% imf为经eemd分解后的各imf分量值
% ICEEMDAN分解
% function imf = pICEEMDAN(data,FsOrT,Nstd,NE,MaxIter)
% 输入：
% data为待分解信号
% FsOrT为采样频率或采样时间向量，如果为采样频率，该变量输入单个值；如果为时间向量，
% 该变量为与y相同长度的一维向量。如果未知采样频率，可设置为1
% Nstd为附加噪声标准差与Y标准差之比
% NE为对信号的平均次数
% MaxIter：最大迭代次数
% 输出：
% imf为经ICEEMDAN分解后的各imf分量值
fs = 100;
t = 1/fs:1/fs:1;
data = sin(2*pi*5*t)+2*sin(2*pi*20*t);
imf = iceemdan(data,0.2,10,100,2);